Regime shifts, trends, and variability of lake productivity at a global scale

Significance Lakes can change dramatically following a slow change in conditions. They can abruptly shift from being oligotrophic to eutrophic or vice versa, in what is called a regime shift. Despite the important consequences for ecosystems and human activities of abrupt shifts, we do not know how frequent they are or how they are distributed globally. To answer these questions, we analyze lake productivity dynamics of 1,015 lakes worldwide. Our results show few experienced regime shifts, yet the occurrence of observed regime shifts is increasing over time. Our analysis' global scope allows us to better understand the occurrence of regime shifts and the socioeconomic drivers associated with them. This knowledge will help manage lakes' response to global change.

individual algorithms' validity ranges. They showed that chlorophyll-a retrieval uncertainties increase below 10 mg/m3, while the detection limit is below 1 mg/m3. We see no other discontinuities than the slightly lower performance at <10 mg. By binning units mg/m3 to TSI, daily to 10-day aggregates, and pixels to lake area averages, we further mitigate the decrease of signal-to-noise for lower concentrations. In other words, every 10-day lake average TSI in our time series corresponds to thousands or millions of pixels, which cancels random noise.

S.1.1.1. HydroLAKES database
General knowledge about lakes, which is crucial for the investigation of ecological processes, was obtained from two different data sources. Our main data source is the HydroLAKES (13), a database of 1.4 million global lakes with a surface area of at least 10 ha, which provides a general description as well as geographic attributes such as location and area of the catchment. It also includes estimates of the lake depth, water volume, and residence time. For lakes where no HydroLAKES information was available, we used the Global Lakes and Wetlands Database (GLWD) (14). This was the case for 15 of 1016 lakes. For 4 lakes, no lake characteristics could be retrieved and were excluded from the statistical analysis.

S.1.3.1. Morphological and geographical characteristics
In order to obtain information at the catchment level, the first step is to delineate the catchments for each of the lakes. The HydroLAKES database mentioned above is a subset of HydroSHEDS, which provides not only the above-described database but also HydroBASINS, a collection of global catchments, crucial for hydrological modeling.
HydroBASINS consists of a variety of layers for the process chain for catchment delineation, from Voidfilled elevation (DEM) to the final basins layer in different resolutions. It represents a seamless global coverage of hierarchically nested sub-basins.
We used HydroBASINS for the catchment delineation of the lakes in our study. Depending on the size or location of the lakes and also considering process time, different approaches are used. We used the flow direction layer in a resolution of 15 arc-seconds (approximately 500m at the equator) as input data for the lake catchment generation. If this layer was not available, we derived the catchment by aggregation of all sub-basins upstream by the hierarchical procedure. Once we obtained the catchment for each lake, we used it as a mask to obtain the catchment-level covariates at the lake level. In particular: human population, precipitation, and sub-national gross domestic product.
All spatial analysis is done in ArcGIS Desktop 10.5 (ESRI 2016) with standard GIS tools. For the catchment delineation, the extension ArcHydroTools in version 10.2 is also used with some minor adaptations in python scripts.

S.1.3.2. Human population
The world population is obtained for the year 2014 from Landscan (15), a global digital raster data in a resolution of approximately 1 km, which is the finest global population distribution data available today.
It was developed for the U.S. Department of Defense and distributed by ORNL's LandScan. We obtain the total number of populations in a catchment as well as the population density (number of people / catchment area in km 2 ).

S.1.3.3. Precipitation and temperature
As data sources for the precipitation, we used the Climatologies at High resolution for the Earth Land Surface Areas -CHELSA in Version 1.1 (16), a global database of the downscaled model output of temperature and precipitation estimates at a horizontal resolution of 30 arc sec (approximately 1 km). The precipitation algorithm includes orographic predictors such as wind fields, valley exposition, and boundary layer height. The resulting monthly precipitation is extracted as time series for a given catchment for each month from 2001 to 2013. The trend in precipitation and the trend in temperature for each catchment are calculated as Spearman's rank correlation between precipitation and time and temperature and time, respectively.

S.1.3.4. Sub-national gross domestic product
The Gross domestic product (GDP) is obtained from the Gridded global datasets for Gross Domestic Product (17). This is a global multiannual dataset, consistent over time and space at 5 arc-min resolution (approximately 10 km) for the 25-year period from 1990 to 2015. The GDP per capita dataset represents the average gross domestic product per capita in a given administrative area unit, whereas, for the total GDP, the above values are multiplied by a gridded population. Both datasets, available as NetCDF files, are retrieved per year and aggregated for each catchment.

S2. Statistical analysis
The variables used in the statistical analysis are reflected in SI Appendix Table S1.

S2.1. Structural equations modeling to determine variables associated with long-term trends
We used structural equations models (SEMs) to test causal hypotheses regarding the effects of climate, geographic, geologic and anthropogenic variables on the mean, variance, and the trends in both the mean and variance of TSI over time. We implemented the SEMs using the pSEM function in the piecewiseSEM package in R v. 4.0.2). We generated causal hypotheses by investigating the variance-covariance matrix of response variables and covariates of interest: mean temperature, mean precipitation, temporal trend in temperature, temporal trend in precipitation, lake depth, catchment area, human density, and the mean GDP of people living within the lake catchment. We began our model selection procedures by testing causal hypotheses for variables explaining the mean TSI. We expected that the mean TSI could also be a predictor for the variance of TSI, the trend in the mean ,and the trend in the variance of TSI. To begin, we used the structured equations to test whether the mean TSI was well-predicted by climate variables (temperature and precipitation), lake properties (lake depth and catchment area), and/or anthropogenic variables (human density and average GDP). Whenever average GDP was an important predictor, we tested whether it, in turn, was well-predicted by average temperature and human density. Similarly, whenever human density was an important predictor, we, in turn, tested whether the average temperature was an important predictor of human density. Average temperature, average precipitation, lake depth, catchment area, human density, and average GDP were all ln(+1)-transformed.
Each equation within an SEM is a linear model. We followed a progression of including missing causal links (identified by tests of directed separation, 'd-sep' tests), as well as removing non-significant pathways, repeating this procedure until all identified significant paths were included. When a path was significant, but we had no hypothesis regarding the causal relationship between the variables, we included the relationship as a correlation. We used Fischer's C statistic and the p-value to evaluate the plausibility of the models. All final models were plausible, with p>0.05.
We started with the following model for the mean TSI:
Starting model for trend in mean TSI:

S3.1. Spatial resolution
There are several limitations to be considered and improved upon in the future. The first one is the spatial resolution of the remote sensing data. This translates into a bias towards large lakes (Fig. S1B) when compared to other reference lake datasets like Hydrolakes. Since lake area and depth are correlated, the spatial resolution limit also translates into a bias towards deep lakes (Fig. S1C), although this bias is not as strong as in the case of lake area.

S3.2. Length of the time series
Even though we are using the longest time series available, they might be too short to show signals of instability in some lakes. For this reason, we cannot describe a lake as stable, only to report the abrupt changes and signs of instability that we observe within the period covered in these time series.

S3.3. Data gaps
For the change-point detection algorithm, a continuous time series sampled at regular intervals is a prerequisite. Meris sensors work on the optical range and cannot penetrate clouds. Therefore, the number of gaps in the time series depends on how much cloud coverage or ice a lake experiences. Lakes in high latitudes are more prone to have long periods without data, particularly during winter. The Spearman's rank correlation between the number of missing data points and the absolute value of latitude is 0.43, with a pvalue < 0.001.
During the gap-filling routine, when the data within a moving window was scarce, the uncertainty around the gap-filling was such that it is increasingly difficult to find significant trends in the data. The correlation between the number of data gaps and the significance of the slope between two change points is -0.47, pvalue < 0.001).

S3.4. Finding signatures compatible with fold bifurcations in observational data
Pure observational data cannot conclusively determine whether a regime shift is the result of a fold bifurcation. Hysteresis could not be detected without knowing what the control parameter might be and quantifying it in situ, the control parameter might not continue to change after the shift, or the system might not relax faster than the sampling frequency. Therefore, we refer to the abrupt changes observed as candidate regime shifts and candidate tipping points.

S3.5. Chosen state variable
Lakes might experience a regime shift in another variable other than TSI within the study period, like turbidity or the population of specific species.    S3. Sensitivity analysis. This figure shows that the frequency of candidate regime shifts from low to high TSI is increasing over time regardless of the number of bins in which we partition the time series. Although not for every partition, the correlation is significant (Spearman's rank correlation p-val < 0.05, orange bars), the pattern is consistent. All correlation coefficients are positive and similar in magnitude. Moreover, this can also be understood as that the time between two candidate regime shifts is decreasing over time. Such correlation (time between candidate regime shifts in the dataset vs. time) would not require a sensitivity analysis. We find that correlation to be significant (p-val < 0.01) Spearman's rank correlation of -0.17, confirming that, over the length of the time series, candidate regime shifts from low to high TSI have become more frequent over time. The temporal trends in the total number of candidate regime shifts (independent of shift direction), and candidate regime shifts from high to low TSI are shown in Fig. S4.  Fig. S4. Panels A and C are equivalent to figure 2C in the main text but for all candidate regime shifts (A) and for lakes that go from high to low TSI as a result of the shift (C). Panels B and D are the sensitivity analysis equivalent to Figure S3. It is clear that the positive (although not significant) trend in the total number of candidate regime shifts is driven by the increase in candidate regime shifts from low to high TSI shown in figure S3. No temporal trends were observed in the number of candidate tipping points.   Table S2: Best model to explain the probability of a lake experiencing a candidate regime shift in TSI (RSTSI is a binary vector indicating whether a lake has experienced a candidate regime shift, 1, or not, 0). After model selection, the final binomial model takes the following form: Table S3: Best model to explain the coefficient of variation in a lake. After model selection, the final model takes the following form: Table S4: Best model to explain the probability of a lake experiencing a candidate tipping point in TSI (TPTSI is a binary vector indicating whether a lake has experienced a candidate tipping point, 1, or not, 0). After model selection, the final binomial model takes the following form: